variables

ts=(t0:stepsize:tend)'; %sampling times
N=length(ts); % Anzahl Schritte

testinp=@(t) sin(omega*t.*t); %input function
u=testinp(ts); %inputs at times ts
disturbance=random('norm',0,ones(size(u))*0.02);

testt=tf([ucoeff;0].',[1;-ycoeff].',stepsize,'Variable','z');
% SYS = TF(NUM, DEN, TS, 'PropertyName1', 'PropertyValue1')
%   creates a discrete-time transfer function with
%   sample time TS (set TS=-1 if the sample time is undetermined), 
%   with numerator(s) NUM and denominator(s) DEN.  
%   The output SYS is a TF object.
%------------------
%testt: 
% 
%Transfer function:
%    0.73 z^3 + 3 z^2 + 7 z
%------------------------------
%z^3 - 0.15 z^2 - 0.27 z + 0.43
% 
%Sampling time: 0.01
%-----------------

y=dlsim(testt,u,disturbance,ts);
% dlsim(sys, u, v, t) returns the time response of the lti model sys', where
% sys is of the form y(z) = B(z)/A(z) * u(z) und sys' is y(z) = B(z)/A(z) +
% 1/A(z) * v(z), that is, a _d_isturbance v (hence dlsim) is added in every
% step. We need this specific model structure because it is the model 
% structure identified by the basic least squares model identification 
% algorithm, that we have programmed and want to test.
% At the moment t is ignored.



% Die Relation y(k) = y(k-1)*a(1)+y(k-2)*a(2)+...+y(k-yl)*a(k-yl) +
% u(k)*b(0)+k(k-1)*b(1)+...+u(k-ul+1)*b(k-ul+1) (+ möglicherweise eine 
% Störung v(k) bestimmt für jedes k eine Gleichung, die linear in [a,b] ist.
% Wir wollen das entsprechende (überbestimmte) lineare Gleichungssystem 
% lösen um [a,b] zu finden.

firstentry=max(yl+1,ul); % Welches k gibt uns die erste lineare Gleichung


% S so zusammengestückelt, dass gilt S*[a.';b.']=ys

ys=y(firstentry:N); 
S = zeros(N-firstentry+1,length(coeff));


for k=1:yl
    S(:,k)=y(firstentry-k:N-k);
end

for k=0:(ul-1)
    S(:,k+1+yl)=u(firstentry-k:N-k);
end


plot(ts(firstentry:N),ys,'m',ts(firstentry:N),u(firstentry:N),'g') %Input und Output des Systems, das wir schätzen

phat=S\ys; 

% Kontrolle: es kommt ca. [ycoeff,ucoeff] raus
phat-coeff
norm(phat-coeff)
